function [y, done] = laguerre ( f, x0, TOL, Nmax )
%LAGUERRE  root refinement method

done = 0;
n = length ( f ) - 1;
fp = polyder ( f );
fp2 = polyder ( fp );

for i = 1 : Nmax
    fx = polyval ( f, x0 );
    fpx = polyval ( fp, x0 );
    fp2x = polyval ( fp2, x0 );

    if ( abs(fx) < TOL ) 
       y = x0; done = i;
       return
    end;

    gx = fpx / fx;
    g2 = gx * gx;
    hx = g2 - fp2x / fx;

    disc = sqrt ( (n-1) * ( n * hx - g2 ) );
    if ( abs(gx-disc) < abs(gx+disc) )
       denom = gx+disc;
    else
       denom = gx-disc;
    end
    dx = n / denom;
    x0 = x0 - dx;

    if ( abs(dx) < TOL ) 
       y = x0; done = i;
       return
    end 
end
y = x0;
end